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■ Abstract. We analyse the evolution of the fractional ionisation in a steady-state protoplanetary disc over 10 6 yr. 

We consider a disc model with a vertical temperature gradient and with gas-grain chemistry including surface 
reactions. The ionisation due to stellar X-rays, stellar and interstellar UV radiation, cosmic rays and radionuclide 
decay is taken into account. Using our reduction schemes as a tool for the analysis, we isolate small sets of chemical 
reactions that reproduce the evolution of the ionisation degree at representative disc locations with an accuracy 
of 50%-100%. On the basis of fractional ionisation, the disc can be divided into three distinct layers. In the dark 
dense midplane the ionisation degree is sustained by cosmic rays and radionuclides only and is very low, S~ 10 -12 . 
This region corresponds to the so-called "dead zone" in terms of the angular momentum transport driven by MHD 
turbulence. The ionisation degree can be accurately reproduced by chemical networks with about 10 species and a 
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similar number of reactions. In the intermediate layer the chemistry of the fractional ionisation is driven mainly by 
the attenuated stellar X-rays and is far more complicated. For the first time, we argue that surface hydrogenation 
of long carbon chains can be of crucial importance for the evolution of the ionisation degree in protoplanetary 
discs. In the intermediate layer reduced networks contain more than a 100 species and hundreds of reactions. 
Finally, in the unshielded low-density surface layer of the disc the chemical life cycle of the ionisation degree 
comprises a restricted set of photoionisation-recombination processes. It is sufficient to keep about 20 species and 
reactions in reduced networks. Furthermore, column densities of key molecules are calculated and compared to 
CZ2 . the results of other recent studies and observational data. The relevance of our results to the MHD modelling of 

protoplanetary discs is discussed. 
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1. Introduction Due to this coupling, the active layer is unstable to the 

. r i i • f ! MRI, and the developing turbulence allows the accretion 

Nowadays a paradigm of the evolution of a protoplanetary , i • w ttm c ch. u\i\iw h 

...... . . ... r i t to occur (e.g., Gammie 1996 Fleming & Stone 2003). 

disc is widely accepted m which most ol the disc matter 



is assumed to move steadily toward a protostar due to re- 
distribution of the angular momentum. Ionisation in such 
objects is an important factor that enables the angular 
momentum transport to occur via magnetohydrodynamic 



The shielded midplane region is almost neutral, decoupled 
from the magnetic field, and, thus, quiescent. 



The location of the boundary between these two re- 

_ rTTT _ , , , . , . gions may prove to be very sensitive to the disc physi- 

MHD turbulence driven by the magnctorotational m- , ,. n , . , , „ 

.... „. T „ ,. „ , ■ , _ , . . cal properties and chemical composition (e.g., hromang, 

stability (MRI; Balbus k Hawley [T» From this pom & Balbus ^ Hq ^ MRB modeUi 

ol view, a disc is conventionally divided into the "active .,, . , , „ ' . n , , . n n . 

J with non-ideal effects included is a very demanding corn- 
layer, adjacent to the disc surface, and the dead zone, , ,. , , , . , . ,, ,, rm , „. c 
, , . „, . putational task. 1ms is why m the MHD modelling of 
centred on the midplane. f he active part of the disc is lr- , ,. / , . . ,, , , >, 

., ii-i.i ii ;• 7 n i mi protoplanetary discs (and protostellar clouds) a very sim- 

radiated by high energy stellar /interstellar photons. Ihus, , , . , , . „ , r 

. . ? . . i . i.ii.i i • i • chemical scheme is usually assumed with a few ions 

the fractional ionisation there is relatively high, which lm- , , . , . . , , , . . . . , 

. ii an d a network that includes only ionisation and recom- 
piles that the magnetic field is well coupled to the gas. , . , . n , . . . , r ^ 
° bmation reactions, often neglecting the presence of dust 

Send offprint requests to: D. Semenov, e-mail: grains (Sanoet al. 120001 Fromang et al. 120021 Fleming 

3ia-hd.mpg.de & Stone l2UH3*|l . The medium is believed to be in chemi- 
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cal equilibrium so that the fractional ionisation x e can be 
expressed as (e.g., Gammie 1996) 



x e (eq) 



C 



(1) 



where £ is the ionisation rate, (3 is a typical recombination 
coefficient, and nn is the hydrogen number density. This 
approach may indeed be valid if only the cosmic ray ion- 
isation is taken into account and only gas-phase chemical 
processes are considered. However, dust plays an impor- 
tant role in the evolution of the fractional ionisation being 
an efficient electron donor for recombining ions and a sink 
for neutrals in cold parts of a disc. Moreover, newly born 
stars possess a relatively high X-ray flux with photon en- 
ergies from about 1 to 5 keV (e.g. Igea & Glassgold 1999). 
These factors give rise to a more complicated chemistry 
relevant for the ionisation degree. In particular, different 
ions dominate the fractional ionisation at different times 
and in different parts of the disc. Therefore, the straight- 
forward application of Eq. Q for evaluating the ionisation 
degree can be fraught with errors in some parts of a disc. 

To check the validity of Eq. JQ), we analyse the de- 
tailed ionisation structure of a protoplanetary disc com- 
puted with the full UMIST 95 chemical network (Millar et 
al. I1997J) with the surface chemistry included. We adopt 
a reference disc model to serve as a guide to the range of 
possible physical conditions that may be encountered in 
real protoplanetary objects. Within this model, we choose 
several representative disc locations and investigate in de- 
tail which chemical processes control the time-dependent 
ionisation degree there. Our intention is to show that the 
oversimplified treatment of the ionisation in a protoplane- 
tary disc as an equilibrated ionisation-recombination cycle 
can lead to x e values that differ by more than an order of 
magnitude from values computed with all the available 
information on the disc chemical evolution. 

We describe the chemical processes that are respon- 
sible for ionisation in a disc in terms of reduced net- 
works that are subsets of the full UMIST 95 network, sup- 
plied in some cases by a few surface reactions. These net- 
works contain only those species and reactions that are 
needed to reproduce x c with up to a factor of 2 uncer- 
tainty. The utilised reduction methods are described in 
Wicbc, Semenov & Hcnning (2003, hereafter Paper I). The 
species-based reduction rests upon the Ruffle et al. ({5002 ) 
technique and consists of choosing species that are im- 
portant in a particular context, and then selecting from 
the entire network only those species that are necessary 
to compute abundances of important species with a rea- 
sonable accuracy. In the reaction-based method, analysis 
starts from reactions that govern the abundance of im- 
portant species. All reactions in the entire network are 
assigned weights according to the influence they have on 
abundances of important species. Then, only those reac- 
tions are selected that have weights above a cut-off param- 
eter that is selected on the basis of the requested accuracy. 



Only those species are included in the reduced network 
that participate in selected reactions. 

The organisation of the paper is the following. In 
Section [2] we describe the disc model and its physical 
structure as well as updates to the chemical model of 
Paper I. In Section we discuss the initial conditions for 
the disc chemistry. The processes responsible for the frac- 
tional ionisation in various parts of the disc are outlined in 
Section^ In Section[S]column densities are tabulated and 
compared to other studies and observational data. Results 
of the analysis and their relevance to the MHD modelling 
are discussed in Section Final conclusions are drawn in 
Section □ 

2. Disc model and physics update 

2.1. Disc structure and parameters 

As a basis for our calculations, we adopted the steady- 
state disc model of D'Alessio et al. ||T599). The central 
star is assumed to be a classical T Tau star with an effec- 
tive temperature = 4 000 K, mass = 0.5M Q , and 
radius i?„ = 2i? Q . The disc has a radius of 373 AU, an ac- 
cretion rate M = 10~ 7 Mq yr _1 , and viscosity parameter 
a = 0.01. It is illuminated by UV radiation from a central 
star and by interstellar UV radiation. The intensity of the 
stellar UV field is described with the standard G factor, 
which is scaled in a way that the unattenuated stellar UV 
field at R = 100 AU is a factor of 10 4 stronger than the 
mean interstellar radiation field (DraineEZHJ- The visual 
extinction of stellar light at a given disc point is calculated 



mag 

1.59- 10 21 cm- 2 ' 



(2) 



where N# is the column density of hydrogen nuclei be- 
tween the point and the star. The extinction of the inter- 
stellar UV radiation is computed in a similar way (but in 
a vertical direction only) . 

In addition to the UV field, three other ionisation 
sources are considered: cosmic ray ionisation (Ccr)> decay 
of radionuclides (Cr), and stellar X-ray ionisation (Cx)- 
The following expression is used to compute the cosmic 
ray ionisation rate: 

Ccr = ( o\ 
±Co [exp (-E x (z, r)/10 2 ) + exp (-£ 2 (*. r)/10 2 )] , W 

where £i(z,r) is the surface density (g cm -2 ) above the 
point with height z at radius r, T,2(z,r) is the surface 
density below the point with height z at radius r, and 
Co = 1-3 • 10~ 17 s _1 . Thus, we assume that cosmic ray 
particles penetrate the disc only in the vertical direction 
from both its surfaces. The X-ray ionisation rate Cx is 
computed according to Glassgold, Najita & Igea (|T997a 
1997b) with parameters for their high depletion case. The 
source of X-rays is assumed to be located at z = 12 Rq. 
We adopt the ionisation rate due to decay of radioisotopes 
to be Cr = 6.5 • 10~ 19 s" 1 . 
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The region under investigation spans a range of radii 
from 1 to 373 AU. We do not consider regions closer to 
the star to avoid the necessity to account for physical pro- 
cesses like dust destruction and collisional ionisation. Even 
when the close vicinity of the star is excluded, the disc is 
still characterised by a wide range of all relevant physical 
parameters: gas temperatures T g are 10 — 600 K, densities 
are 10~ 20 - lO" 10 g cm" 3 , G are 10 2 - 10 8 , Ay varies 
from to more than 100 mag and ionisation rates are 
10-18 _ iQ-io s-i. 

We divide the disc into three regions, namely, the mid- 
plane, the intermediate layer, and the surface layer. This 
division is similar to that outlined by Aikawa & Herbst 
(1999}. The midplane is the part of the disc near the 
symmetry plane, which is opaque to both X-ray and UV 
photons and is ionised primarily by cosmic ray particles. 
If one goes up from the midplane, the fractional ionisa- 
tion first stays the same as in the midplane, but then at 
some height it starts to grow in response to increasing X- 
ray intensity. We take this turn-off height to be the lower 
boundary of the intermediate layer (and the upper bound- 
ary of the midplane region). In our model this height is 
approximately given by zl = 0.05 r 142 AU. With increas- 
ing height, the fractional ionisation grows first slowly, then 
much faster, in response to the decreasing opacity of the 
medium. The point where transition between slow and 
fast growth occurs is assumed to be the upper boundary 
of the intermediate layer and the lower boundary for the 
third, surface layer. This boundary is approximated by 
z\j = 0.17r 127 AU. The fractional ionisation is important 
in dynamical calculations only if it is small, as high ioni- 
sation means that the ideal MHD equations can be used. 
Thus, we set the upper boundary of the surface layer to 
a position where the ionisation degree exceeds ~ 10~ 4 . 
This simple picture may not be appropriate if diffusion 
processes or radial transport of the disc matter are taken 
into account (Ilgner et al. 2003). 

The disc structure is depicted in Figure ^ The up- 
per boundary of the surface layer is shown schematically 
(thick solid line), as at R 100 AU it extends beyond the 
disc limits (dotted line) for the adopted cut-off fractional 
ionisation of 10~ 4 . In Figure[21we show the dependence of 
the fractional ionisation on the height above the disc plane 
computed with the chemical model, described below, for 
several representative radii R. 

2.2. Chemical model 

The chemical model used in this paper is essentially the 
same as in Paper I, but with a few corrections that are de- 
scribed in this subsection. We adopt the UMIST 95 rate- 
file (Millar et al. I1997J) for the gas-phase chemistry and 
the surface reaction set compiled by Hasegawa, Herbst 
& Leung i 1 99211 and Hasegawa & Herbst (1993;). As we 
have multiple ionisation sources in the protoplanetary 
disc, the cosmic ray ionisation rate is replaced by the sum 
Ccr + Cx + Cr m expressions for rates of chemical reactions 




100 200 300 400 

R, AU 

Fig. 1. Three layers of a disc with distinctively different 
sets of chemical processes responsible for the value of the 
fractional ionisation x c . The thick solid line is the upper 
boundary of the surface layer where x c = 10~ 4 , while 
dotted line depicts the disc upper limits. 
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Fig. 2. Fractional ionisation as a function of height above 
the disc plane. Lines are labelled with the corresponding 
radii. 

with cosmic ray particles (CRP) and cosmic ray induced 
photoprocesses. Given the high intensity of the impinging 
UV radiation, photodesorption of surface species is taken 
into account along with other desorption mechanisms. 

2.2.1. Grain charge 

The major difficulty in modelling the fractional ionisation 
in a dark and dense environment, like a protoplanetary 
disc midplane, is that grains play a very important role 
in ion recombination and can be the dominant charged 
particles in some disc regions. The charge of grains and, 
consequently, their impact on the local dynamics can sen- 
sitively depend on the grain size distribution, icy mantle 
composition etc. (e.g. Nishi, Nakano & Umebayashi 1991). 
Thus, the computation of the fractional ionisation in the 
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disc midplanc is related not only to the chemistry but to 
the grain physics as well. This means we should include 
at least a simplified treatment of grain charging processes 
in the model. 

The following changes are made with respect to 
Paper I. We consider neutral, positively charged, and neg- 
atively charged grains separately. A grain attains a neg- 
ative charge or loses a positive charge colliding with free 
electrons. The probability of an electron sticking to grain 
surfaces is assumed to decrease exponentially with dust 
temperature Td, so that the sticking coefficient is around 
0.5 at Td ~ 10 — 20 K, close to the value obtained by 
Umebayashi & Nakano (TX980), and is essentially zero at 
higher temperatures. Via dissociative recombination of a 
gas-phase ion, a grain loses one electron and becomes neu- 
tral or positively charged. Rates of collisions between elec- 
trons and positively charged grains are multiplied by a 
factor, correcting for the long-distant Coulomb attraction, 
Cion = 1 + e 2 /fcadTd, where ad = 0.1 (im is the grain ra- 
dius. It is always assumed that T g = Td = T '. 

2.2.2. Sticking probability 

In Paper I the sticking probability S was taken to be 0.3 
at T g = 10 K for all neutral species except H, He, and H2. 
However, at higher temperatures S is likely to be smaller 
(Burke & Hollenbach 1983). To account for this tendency, 
we multiply S for neutrals by an additional factor equal 
to the fraction of molecules of a given kind that have ther- 
mal energy lower than the desorption energy Ed for this 
species (assuming a Maxwellian velocity distribution). 

2.2.3. Desorption processes 

The disc model, used in this paper, is characterised by a 
much wider range of physical conditions than the model 
for a molecular cloud considered in Paper I. This means, 
in particular, that changes must be made to the adopted 
model of the species desorption from grain surfaces. First, 
we now have the UV radiation in the model, so we add 
photodesorption to the processes listed in Paper I. The 
rate of photodesorption is given by 

fc ph = [G s exp(-2^) + Gis ex p(-2^ s )] Ynal (4) 

where Gs and Gis are intensities of the stellar and inter- 
stellar UV radiation expressed in units of the mean inter- 
stellar radiation field, A v and ^4y are the corresponding 
visual extinctions in the direction to the star and in the 
vertical direction, and Y is the photodesorption yield for 
which we adopted the expression 

Y = 0.0035 + 0.13 exp(-336/T d ), (5) 

derived by Walmslcy, Pincau des Forets & Flower (1999) 
from the experimental data obtained by Westley et 
al. (JTM3- 

The cosmic ray desorption rate also needs to be 
revised. The expression suggested by Hasegawa & 



Hcrbst (1993J and used in Paper I is based on the as- 
sumption that a cosmic ray particle (usually an iron nu- 
cleus) deposits on average 0.4 MeV into a dust grain of the 
adopted radius, impulsively heating it to a peak temper- 
ature T crp , which is close to 70 K for cold 0.1/mi silicate 
dust. As in a disc we expect much higher dust tempera- 
tures, data from Leger, Jura & Omont l|1985fl are used to 
take into account the dependence of the cosmic ray heat- 
ing on the initial dust temperature. The peak temperature 
is approximated as 

T crp = (4.36-10 5 + T d 3 ) 1 / 3 . (6) 

This expression predicts that a grain heats up to 76 K at 
Ti = 10 K and gives almost no heating for Td 100 K. 

3. Initial conditions for chemistry 

The problem of the initial conditions for chemical models 
of star-forming regions remains an open issue. The mod- 
elling of pre-protostellar objects commonly starts with the 
purely atomic or partly ionised initial composition, while 
in reality the initial abundances even for a moderately 
dense medium must reflect its previous evolution which 
ideally would have started from a diffuse intercloud gas. 
Even though the atomic or ionic initial composition can 
be adequate in the prestellar or pre-protostellar objects, 
the situation is different in protoplanetary discs, where the 
chemical composition at the beginning of the modelling is 
obviously far more advanced. One way to cope with this 
problem is to use "sliding" initial conditions, when the ini- 
tial abundances are set to be mostly atomic at the outer 
radius of the disc. Then, a chemical model is run, and the 
final abundances at a given radius are used as initial abun- 
dances at the next radius, closer to the disc centre (Bauer 
et al. 119971 Willacy et al. I1998|) . Another approach is to 
compute the input abundances with the model of a dark 
cloud, out of which the disc has evolved (e.g. Aikawa & 
Herbst [TUM Aikawa et al. IMrZj) . Willacy et al. ifTMSjl 
compared both methods and found that, with a few ex- 
ceptions, they provide similar results, because in a high- 
density environment reasonable molecular abundances are 
reached within a fraction of a year even with the atomic 
initial composition. 

To simplify the reduction, we adopted the second ap- 
proach. The choice of a cloud model in Aikawa et al. ( 2002 ) 
and in previous papers of this group is based on the 
analysis of the Ohio New Standard Model performed by 
Terzieva & Herbst (|1998al I1998b|l . These authors devel- 
oped a simple composite criterion that allows one to es- 
timate how good the chemical model is in reproducing 
the chemical composition of a given object. Adopting a 
fixed density of 2 ■ 10 4 cm~ 3 and a gas temperature of 
10 K, Terzieva & Herbst l|1998a|> found that with the pure 
gas-phase chemistry, the chemical composition of the two 
typical molecular clouds, namely, TMC-1 and L134N, is 
best reproduced after about 3 • 10 5 years of the evolution. 
Aiming to reproduce the observed low abundances of wa- 
ter and molecular oxygen, Roberts & Herbst H2002|l added 
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Table 2. Physical conditions in the midplanc. 



3 4 5 6 
Log Time, Yrs 



Fig. 3. Percentage agreement between calculated and ob- 
served gas-phase abundances. 



surface chemistry to this model and found that the "best 
consistency" interval can be shifted toward a somewhat 
later time, like a few times 10 6 years. 

As we use the UMIST95 ratefile, we perform a sim- 
ilar analysis for this database. The observed abundances 
for TMC-1 and L134N are taken from Ohishi, Irvine & 
Kaifu Ijl992[l and Langer et al. (1997). Following Terzieva 
& Herbst f 199S1J. we assume that the abundance of a 
given molecule is consistent with observational data if it 
lies within an order of magnitude of the observed abun- 
dance. When the observed abundance is an upper limit, 
we assume that the computed abundance is in agreement 
with observations if it is less than or equal to 10 times this 
upper limit. With the same density, temperature and the 
surface chemistry included, we find that for our chemical 
model the agreement is best at t ~ 10 6 years (Fig. [31) ■ The 
"low-metal" set of abundances from Paper I is adopted for 
this computation. 

We now use the set of gas-phase and surface abun- 
dances, obtained with the above model for t = 10 6 years 
as the initial conditions. The brief outline of this set is the 
following (abundances relative to the number of hydro- 
gen nuclei are given in parentheses): The primary carbon- 
bearing compound is surface formaldehyde (3.7 • 10~ 5 ) 
which also accounts for a significant fraction of the over- 
all oxygen content, being second only to surface water 
(7.8T0 -5 ). Gas-phase oxygen and carbon are locked in CO 
molecules (2.8-10 -5 ). A fraction of oxygen is also present 
in atomic (1.1 • 10 -5 ) and molecular (8.7 • 10~ 6 ) form. 
Gas-phase nitrogen is contained mainly in N2 molecules 
(9.7 ■ 10 -6 ), while in icy mantles the main N-bearing 
species are HCN and HNC (~ 10 -6 each) with a slightly 
lower amount of CH2NH, NH3, and NO. Sulphur is dis- 
tributed almost equally between gas and dust phases. In 
the gas-phase it is locked in the CS molecules (2.9 • 10~ 8 ). 
Surface sulphur is bound in H 2 S (2.2 • 1CT 8 ) and H 2 CS 
(1.5 • 10~ 8 ). Even though the density is not very high, 
metals are already significantly depleted. Surface abun- 



No. 


R, AU 


p, g cm :i 


T, K 


CCR, 


s- 1 


Cx, 


s" 1 


Ml 


1. 


2.1(-10) 


614. 


8.8(- 


-19) 


7.5(- 


-18) 


M2 


3. 


4.7(-ll) 


193. 


2.3(- 


-18) 


4.0(- 


-23) 


M3 


10. 


8.5(-12) 


52. 


5.0(- 


-18) 


1.2(- 


-25) 


M4 


30. 


7.3(-13) 


29. 


9.4(- 


-18) 


5.4(- 


-27) 


M5 


100. 


4.8(-14) 


16. 


1.2(- 


-17) 


3.8(- 


-28) 


M6 


300. 


4.8(-15) 


9. 


1.3(- 


-17) 


4.0(- 


-29) 


M7 


370. 


3.3(-15) 


8. 


1.3(- 


-17) 


2.6(- 


-29) 



Table 3. Reduced network for dark, hot chemistry in the 
midplanc. 



NO + + Mg 
Mg + + e~ 
NO + hven 



Mg + + NO 
Mg + hi; 
NO+ + e" 



dances of Mg, Fe, and Na constitute about 1/3 of their to- 
tal abundances. The dominant ions are HCO + (3.8 • 10~ 9 ), 
H+ (3.5 • 1(T 9 ), Mg+ (1.9 • 1(T 9 ), C+ (1.3 • 1CT 9 ), and Fe+ 
(1.2- 10- 9 ). 

It should be kept in mind that some of our reduced 
networks computed for this initial composition may not 
be appropriate for purely atomic initial abundances of for 
abundance sets with higher metal content. 

4. Chemistry of the ionisation degree in a 
protoplanetary disc 

In each of the selected disc points the chemical model 
described in Section F2] is run up to t = 10 6 years. 
Dominant ions in various disc domains are listed in 
Table ^ Reduction methods are then used to find those 
processes that determine x c . All the reduced networks de- 
scribed in this section are freely available from the au- 
thors. Also, we show some representative networks for the 
midplane and the surface layer. 

4.1. Ionisation in the midplane 

The midplane of the disc is characterised by high density, 
high optical depth and a relatively low ionisation rate. 
Neither the UV radiation from a central star nor the inter- 
stellar UV radiation is able to penetrate into this region. 
The ionisation rate is close to 10 -17 s _1 . At the point 
closest to the star, £ is dominated by X-ray ionisation; 
further away from the centre the attenuated cosmic rays 
are the only ionising factor. The temperature drops from 
~ 600 K at R = 1 AU to about 8 K at the disc edge. 
In the same distance interval, the gas density goes down 
from 10~ 10 g cm~ 3 to less than 10~ 14 g cm -3 . Detailed 
physical conditions for the midplane are given in Table |21 



4.1.1. Dark, hot chemistry 

The chemistry of the point Ml (see Table EJ) is extremely 
simple. In essence, it involves only five species, which are 
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Table 1. Dominant ions in the midplanc, intermediate layer and surface layer at t = 10 6 years. 



R, AU 


1 


3 


10 


30 


100 


300 


373 


Midplane 


Na+ 


HCNH+ 


HCO+ 


HCO+ 


N 2 H+ 


H 3 + 


H 3 + 


Intermediate 


Mg+ 


HCO+ 


HCO+ 


HCO+ 


HCO+ 


HCO+ 


HCO+ 


layer 




S+ 
































and others 












Surface 


C+ 


C+ 


c+ 


C+ 


C+ 


C+ 


C+ 


layer 


H+ 


H+ 













neutral and ionised magnesium, neutral and ionised NO, 
and electrons. The reduced network for this point is shown 
in Table |3J At the beginning of the computation the elec- 
tron abundance decreases from the initial value of ~ 10 -9 
to the new equilibrium value of about 10 -11 , determined 
by the electron exchange between Mg, Mg + , NO, and 
NO + , and then remains nearly constant. With the full 
network, the dominant ion is Na + , but it serves mainly as 
a passive sink for positive charges. This is why it is taken 
out of the reduced network. It suffices to include ionised 
magnesium, which is one of the most abundant ions in the 
initial abundance set. 

Note the apparent similarity of this network to 
the one that determines x e in dense interstellar clouds 
(Oppenheimer & Dalgarno ll974|l in the sense that it con- 
sists of a representative metal, a representative molecu- 
lar ion and the charge transfer reaction between them. 
The explanation is that, even though we include gas- 
grain interaction in our model, at such high temperature 
the chemistry mostly proceeds in the gas-phase. Thus, 
the ionisation degree is determined by gas-phase reac- 
tions. However, due to much higher density these reac- 
tions arc different from those described in Oppenheimer 
& Dalgarno l|1974[) . The major difference is that the re- 
action of H2 ionisation by cosmic rays is not important 
in the considered point. A modest electron supply is pro- 
vided by NO ionisation by cosmic ray induced photons. 
The importance of the NO + ion is determined by its chem- 
ical inertness. Due to this inertness, most of the time its 
abundance exceeds that of other ions, involved in charge 
transfer reactions with magnesium, by more than an order 
of magnitude. With this 5x3 network, the error in the 
fractional ionisation does not exceed 40% during the entire 
computation time. By adding a few more species to the 
reduced network, it is possible to decrease this uncertainty 
to 20%. These species account for the rapid destruction of 
the two other initially abundant ions, i.e. HCO + and H^~. 

The chemistry of the point M2 is also extremely simple, 
but in a totally different way. At this somewhat lower tem- 
perature, T ~ 200 K, magnesium is irreversibly depleted 
onto dust grains. While at earliest times (t < 3 years) 
the evolution of the electron abundance is governed by 
recombination of Mg + , later the fractional ionisation is 
equilibrated and determined by a simple network involving 



H, H2, , and HJ . Included reactions are H2 ionisation 
by cosmic rays, H^" formation, Hjj" and Mg + recombina- 
tion, Mg accretion onto dust grains, and grain (re)charging 
processes. The error in the fractional ionisation does not 
exceed 10% during almost the entire computation time, 
reaching 40% only at the very beginning of the computa- 
tion due to neglected Na + and Fe + ions. Note that the 
equilibrated gas-phase electron abundance x c ~ 10 -13 at 
this point is several times lower than the abundance of 
charged grains. Thus the net charge density at these con- 
ditions is determined not by the chemical kinetics, but 
mainly by grain charging processes. 



Dominant ions in a model of a disc interior which is 
similar to the one used here are determined by Markwick 
et al. (|2002). In the second part of their Table 1 they 
present ions that define the x e value for various heights 
and at R < 10 AU, when X-rays are taken into account. 
The difference to our results is obvious. In our model, 
the dominant ions throughout the inner part of the mid- 
plane are Na + , HCNH+, and HCO + , while in the model 
by Markwick et al. (|2002|l everywhere in the midplane the 
fractional ionisation is determined by the rather complex 
CH 3 CO + ion. We attribute this discrepancy to the differ- 
ent treatment of the surface dissociative recombination of 
ions as well as to different treatment of grain charges. In 
the Markwick et al. model all grains are assumed to par- 
ticipate in recombination reactions, with products of re- 
combination sticking to grain surfaces, while in our model 
these products are returned to the gas immediately. So, in 
their model, metal ions probably recombine with grains 
and stick to dust surfaces, which implies that they are not 
easily returned to the gas-phase. In our model, in the point 
closest to the star, all dust grains are positively charged, 
so they do not contribute to the recombination rate at 
all. Also, when a metal ion happens to encounter a rare 
negative or neutral grain, it loses the charge but remains 
in the gas-phase. Higher above the midplane and some- 
what further from the star, where details of the gas-grain 
interaction are less important, the dominant ions in both 
models are almost the same, i.e., HCNH+ and HCO + . 
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Table 4. Reduced network for dark, cold chemistry in the 
midplane. 



Ho + CRP 


_> jj + _|_ e - 


+ H 2 


-»■ H+ + H 


H+ 4- No 


_ > NoH + 4- Ho 


H+ + e~ 


H 2 +H 


HCO+ + e" 


-> CO + H 


N 2 H+ + e~ 


-> N 2 +H 


Fe+ + g~ 


-> Fe 


Mg+ + g- 


Mg 


h+ + <t 


-> H 2 +H 


HCO+ + g- 


CO + H 


N 2 H+ + g~ 


-» N 2 + H 



4.1.2. Dark, cold chemistry 

The ionisation at the points M3 and M4 displays sim- 
ilar behaviour. In a cold and dense environment, metals 
are depleted onto dust grains almost immediately, and the 
gas-phase fractional ionisation drops to a very low equilib- 
rium value - about 3 • 10~ 15 in M3 and about 8 • 10~ 14 in 
M4. These values are determined by three almost equally 
important cycles, ionisation and recombination of helium, 
formation and destruction of HCO + and formation and 
destruction of Hjj\ The reduced network, consisting of 
these species and processes, reproduces the fractional ion- 
isation with a few per cent accuracy in all but the earli- 
est times. The electron abundance at t < 0.1 yrs in the 
absence of metals is determined by a more complicated 
set of processes, with the fractional ionisation being re- 
produced within 50% uncertainty by a network of about 
40 species. Again, the total charge density at these points 
is determined by charged grains, not by gas-phase ions 
and electrons. 

The dominant ions at points M5, M6, and M7 are 
N 2 H+, HCO+ and H^. The reduced network for the en- 
tire outer part of the disc consists of about 25 species and 
of a comparable number of reactions, which govern abun- 
dances of these ions. At earlier times Mg + and Fe + are im- 
portant as electron suppliers, so in addition to N 2 H + and 
HCO + chemistry, the reduced network includes also neu- 
tral and ionised metals. Among the included reactions are 
metal sticking to grains, cosmic ray ionisation of molecular 
hydrogen, formation, charge transfer between metals 
and H J , and dissociative recombination of ions with grains 
and electrons. The network responsible for a; e at the outer 
disc edge is presented in Table^J It also includes reactions 
of adsorption and desorption of neutral components that 
are not shown for brevity. 

4.2. Ionisation in the intermediate layer 

The intermediate layer differs from the midplane by a 
much higher X-ray ionisation rate. Ranges of physical 
conditions for the intermediate layer are summarised in 
Table They imply that many molecules are abundant 
in the gas phase in this disc domain due to shielding of the 



Table 5. Physical conditions in the intermediate layer. 
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10 
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-17) 




-15) 
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-17) 


-T(- 


-17) 
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1(- 
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-1(- 


-17) 


20 - 40 


1(- 


-17) 


_4(- 
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UV radiation and effective desorption by X-rays. From the 
chemical point of view, the characteristic feature of the in- 
termediate layer is that a chemo-ionisation equilibrium is 
sometimes not reached during the entire 1 Myr time span 
considered. 

4.2.1. Warm, X-ray driven chemistry 

Similar to the midplane, the part of the intermediate layer 
closest to the star (R < 3 AU) is characterised by tem- 
peratures of a few hundred K. As the grain icy mantles 
are nearly absent, the chemistry in this region is driven 
by gas-phase reactions. Dominant ions are metals, HCO + 
and HCNH + . Apart from these components, the reduced 
network contains CO and neutral metals as well as neu- 
tral and ionised N and neutral O along with their hydrides 
[N,0]H„ and [N,0]H+, involved in synthesis and destruc- 
tion of the hydrogen-saturated molecules NH3 and H2O. 
Also included are N 2 , He, He + , which are needed closer to 
the bottom of the layer (where T is high and £ is relatively 
low), and HCN, HNC, HCN+, 2 , and O+ which are im- 
portant further away from the midplane (where T is low 
and £ is very high). The fractional ionisation computed 
with this reduced network differs from the one computed 
with the full network by less than 40% during the entire 
time span. 

Again comparing our data for the intermediate layer 
to those by Markwick et al. (2002) we find that at 3 < 
R < 10 AU our results generally agree, especially given 
the fact that the sampling of the inner part of the disc is 
much more detailed in the Markwick et al. model. Again, 
in the innermost point we find a dominant metal ion in our 
model compared to molecular ions HCNH + and H4C 2 N + 
in theirs. We believe that the origin of this difference is 
the same as described above. 

4.2.2. Cold, X-ray driven chemistry 

The ionisation degree further out from the star (at R > 
30 AU) is established by the evolutionary sequence which 
consists of the following stages. Initially, the more abun- 
dant HCO + and species are destroyed almost imme- 
diately and attain abundances of the order of 10~ n at 
t ~ 10~ 3 years. Metal ions are neutralised somewhat more 
slowly through recombination with negatively charged 
grains and electrons and stick to dust surfaces. The overall 
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electron concentration a; c decreases from ~ 3 • 10 9 cm 3 
at t = to 10~ 9 cm~ 3 at t ~ 0.3 yrs. 

Then, few equilibrium states are reached with a se- 
quence of molecular ions, mainly, NHj, H 3 CO + , and 
HCO + . Equilibrium abundance of the next ions in this 
sequence is higher and is reached later, thus, the time 
dependence of the fractional ionisation has a step-like ap- 
pearance (after the initial decrease). The "second birth" 
of HCO + is caused by H2CO coming from icy mantles. 
Formaldehyde reacts with to form H3CO + , which is 
the dominant ion at 0.1 < t < 30 yrs. One of the recom- 
bination channels for H3CO + maintains its equilibrium, 
restoring H2CO. The other two channels lead directly and 
indirectly to CO production and eventually increase the 
abundance of this molecule up to the point where a new 
equilibrium state is reached for HCO + . This ion eventu- 
ally dominates the fractional ionisation in the outer part 
of the intermediate layer. Direct formation of HCO + in 
the reaction of Hjj~ with CO is damped initially in favour 
of the reaction with H2CO both because of the lower 
initial CO abundance and because the latter reaction has 
a higher rate coefficient a — 6.3 • 10~ 9 cm 3 s _1 compared 
with a = 1.7 • 10~ 9 for H+ + CO reaction. 

To reproduce the step- like x c behaviour, one needs a 
reduced network that includes a few tens of species in- 
volved in cycles of synthesis and destruction of the above 
ions. This network reproduces the nature of the equilib- 
rium stages and their length with an uncertainty that does 
not exceed 20%. 

This simple picture is not appropriate closer to the star 
(3 < R < 30 AU), where the density is high enough and 
the temperature is low enough to allow effective gas-dust 
interaction. At the same time, a high ionisation rate leads 
to an increased abundance of helium ions, which signifi- 
cantly alters the late-time evolution of the fractional ioni- 
sation. An example of this evolution is shown in Figure 0] 
It corresponds to the case with rtn = 8.4 • 10 10 cm~ 3 , 
T = 74 K, C = 2 ■ 10~ 13 . The "exact" solution is shown 
with the solid line. The initial decrease of the electron 
abundance is followed by a few steps corresponding to dif- 
ferent equilibrium states. The last one, that of HCO + , is 
reached after approximately 300 yrs of evolution. If there 
was no relatively abundant ionised helium, the fractional 
ionisation would stay constant from that moment. It is 
this ion, coupled to high density and low temperature, 
that causes a new development. At t ~ 3 • 10 4 yrs the 
ionisation degree increases by about an order of magni- 
tude up to ~ 10~ 8 cm~ 3 , and another equilibrium state 
is reached, with S + being a dominant ion and a dominant 
sulphur-bearing species. 

In the initial abundance set sulphur is mostly bound in 
CS and H2S molecules, which are destroyed by He + , pro- 
ducing either S + or S (as well as C, C + , and H2). Passing 
through a few reactions, the most important of which is 
the reaction with CH4, sulphur is reunited with carbon in 
a CS molecule. This process is well equilibrated, when the 
abundance of methane is high, and leads to a negligible 
S + abundance at t < 10 4 yrs. However, with time, some 
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Fig. 4. Evolution of the fractional ionisation in the inter- 
mediate layer at R = 3 AU. The exact solution is shown by 
the solid line. The dashed line corresponds to the reduced 
network consisting of 112 species and 195 reactions. The 
dotted line (right Y-axis) shows the relative uncertainty in 
the fractional ionisation computed with the reduced net- 
work. 



carbon atoms (both free and incorporated into CH4) are 
taken out of this cycle due to formation of long carbon 
chains in reactions of C + with methane and other hydro- 
carbons. After 3000 yrs the cycle is broken entirely. As 
the abundance of methane and other light carbon-bearing 
species goes down, that of ionised sulphur grows, causing 
almost an order of magnitude increase in the fractional 
ionisation. 

Further evolution of x c is implicitly controlled by gas- 
phase and surface reactions involving carbon chains. The 
reduction methods allow one to find which of them are 
most important, as shown in Figure [S] The main route 
for synthesis and destruction of long carbon chains in the 
depicted scheme starts with an acetylene molecule con- 
taining only two carbon atoms, which is rapidly converted 
to C5H2 molecules. This molecule then either transforms 
to even longer chains, or proceeds further down the cy- 
cle to CsII + which recombines to C4H. This species starts 
a chain of very slow reactions with atomic oxygen that 
degrades it to C3H, C2H, and CH. Eventually, destruc- 
tion of carbon chains leads to the new growth of methane 
abundance and to the slow decrease of the fractional ion- 
isation. Ionised sulphur is consumed in the reaction with 
CH4 and eventually is converted back to CS. The abun- 
dance of H^" determines the ionisation degree at the end 
of the computation. 

Conversion of C5H2 molecules into longer chains (lower 
part of Figure EJ) slows down the methane abundance 
restoration. The reduced network must include almost 
all species with six carbon atoms to reproduce properly 
the time scale of this process. Three key reactions of this 
carbon chain cycle proceed on dust grain surfaces. These 
are reactions of hydrogen addition to C5H, C6, and CeH 
(shown with dashed lines). Without these reactions the se- 
quence of carbon chain synthesis would be terminated by 



D. Semenov et al.: Reduction of chemical networks. II 



9 



C N 

5 



2 2 



CH ►CH 

4 2 5 
I 



C 2 H 



Table 7. Reduced network for X-ray dominated chemistry 
in the surface layer (superscript "d" is used to denote sur- 
face species). 




6 2 



Fig. 5. Main routes governing synthesis and destruction 
of long carbon chains. Surface reactions are indicated by 
dashed lines. 



the C5N molecule that lacks effective destruction path- 
ways in the given conditions. Exhaustion of atomic and 
ionised carbon would keep the ionised sulphur at a high 
level, reached at t ~ 3000 years. Thus, we must conclude 
that at least within the limits of the adopted chemical 
model some surface reactions are of crucial significance. 
As far as we know, this is the first indication that the 
surface chemistry may have an "order-of-magnitude" im- 
portance for the evolution of fractional ionisation in inter- 
stellar matter. 



4.3. Ionisation in the surface layer 

The surface layer encompasses those disc regions where 
the degree of ionisation reaches a value of ~ 10~ 4 at the 
end of the evolutionary time. The corresponding physi- 
cal conditions imply low densities (p ~ 10~ 19 — 10~ 13 g 
cm~ 3 , moderate temperatures (T < 250 K), a high X-ray 
ionisation rate and low obscuration of the UV radiation. 
Conditions in the surface layer are given in Table Note 
that the temperature of the disc atmosphere is nearly con- 
stant in the vertical direction at any distance from the 
star. This rather extreme environment presumes a chemi- 
cal evolution typical of photon-dominated regions (PDR) . 
The chemistry is dominated by gas-phase processes, while 
the role of gas-grain interactions is negligible, apart from 
the desorption of initially abundant surface species. As our 
definition of the surface layer is based on the fractional 
ionisation value, the layer can be further divided into two 
chemically distinct zones. Closer to the star (r SS 10 AU) 
x e is controlled by the X-ray ionisation, while in the outer 
regions (r <; 1Q AU) ionisation is mainly provided by the 
UV photons. 
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4.3.1. X-ray dominated chemistry 

The evolution of the ionisation degree in the hot (T ~ 
200 K) part of the surface layer closer to the star is de- 
termined by X-ray ionisation of atomic hydrogen. The fi- 
nal equilibrium ionisation degree is well reproduced by 
the network, consisting of hydrogen-bearing species (H, 
H + , H^, H3 ) only. This equilibrium value is reached after 
100 years of evolution. Up to this time x e is determined 
by a more complicated set of chemical processes, which in- 
cludes H^" formation, an H 3 + formation and destruction 
cycle as well as reactions involving CO and HCO + . The 
primary carbon and oxygen carrier, formaldehyde, adds to 
this complexity by being desorbed from grain surfaces and 
destroyed by cosmic ray and X-ray induced photons. The 
cycles that involve these molecules are not equilibrated. 
Most reactions with polyatomic species tend to produce 
atomic hydrogen, which is gradually ionised by intense 
X-ray radiation. Eventually, all complex molecules are de- 
stroyed, and at the end of the computation, the gas mostly 
consists of atoms and atomic ions. The entire reduced net- 
work comprises 21 species and 27 reactions and is pre- 
sented in Table [7| The uncertainty does not exceed 60% 
during the entire computational time (Fig. |HJ). Without 
this simple, but important ion-molecule chemistry, the un- 
certainty at earlier times would be much larger (exceeding 
an order of magnitude). 
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Table 6. Physical conditions in the surface layer. 
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Table 8. Reduced network for UV-dominated chemistry 
in the surface layer. 
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Fig. 6. Evolution of the fractional ionisation in the surface 
layer at R — 1 AU. The exact solution is shown with the 
solid line. Dashed line corresponds to the reduced network 
consisting of 21 species and 27 reactions. Dotted line (right 
Y-axis) shows the relative uncertainty in the fractional 
ionisation computed with the reduced network. 

4.3.2. UV-dominated chemistry 

The evolution of the ionisation fraction in the distant re- 
gions of the surface layer (r <; 10 AU) is similar but more 
complicated from the chemical point of view than in the 
case of the X-ray dominated chemistry. As the intensity 
of stellar X-rays decreases, UV photons become the main 
ionising source (see Table . The typical gas tempera- 
ture in this part of the disc decreases below 100 K. The 
entire ionisation evolution of the region is determined by 
formaldehyde desorption and destruction. Icy mantles illu- 
minated by UV radiation evaporate, thus delivering H2CO 
into the gas phase. This molecule is either destroyed di- 
rectly by photons (H2CO — > CO,HCO + ) or ionised and 
dissociated (H 2 CO -> H 2 CO+ -> CO). HCO+ dissocia- 
tively recombines to CO as well, and this molecule is dis- 
sociated by photons producing C and O. The destruction 
processes proceed very fast, so after only a fraction of a 
year the gas composition is mostly atomic, with C + being 
the dominant ion. When all polyatomic species are de- 
stroyed, the fractional ionisation is regulated by an equi- 
librium ionisation-recombination cycle of atomic carbon. 
The reduced network consisting of less than 15 species and 
about 20 reactions leads to less than 20% uncertainty dur- 



ing the entire computational time. An example of such a 
network for R — 100 AU is given in Table [8] 

5. Column densities 

A useful product of our study is the complete chemical 
structure of the disc. We present the calculated column 
densities of some species along with observed values and 
column densities obtained by other groups in Table [9j 
The observed values are taken from Aikawa et al. ( 2002 ) , 
Table 1. Five theoretical models of the chemical evolution 
of a protoplanetary disc are considered. The main differ- 
ences between them are discussed below. 

Willacy & Langer lj2"UUU|> used the gas-phase UMIST 95 
database, supplied with an additional set of gas-grain 
processes and surface reactions. They adopted an ex- 
trapolated flared T Tau disc model of Chiang & 
Goldreich (1997), taking into account cosmic rays and 
stellar UV radiation and assuming a sticking probabil- 
ity S = 1. In all the other models, namely, in Aikawa 
& Herbst flMfl ISnUTf t. Aikawa et al. lj2UU2jl . and van 
Zadelhoff et al. PUT?3|) . the New Standard Model (NSM) 
with a set of gas-grain interaction processes is used as the 
chemical network. 

Aikawa & Herbst l|19991 |2001) considered an 
extrapolated minimum-mass solar nebula model of 
Hayashi (1981), taking into account stellar X-ray radia- 
tion in addition to the UV radiation and cosmic rays and 
assuming an artificially low sticking probability S = 0.03 
that compensates for the absence of non-thermal desorp- 
tion mechanisms. The flared T Tau disc model of D'Alessio 
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et al. I|1999|) for various accretion rates is used with a stick- 
ing probability S — 1.0 in Aikawa et al. (|2002ll and van 
Zadelhoff et al. I|2003[l . Stellar X-rays are not included in 
the latter two models. 

Column densities of selected species presented in these 
papers are compared to those computed in this work in 
Table Wherever possible, we took values from tables, 
otherwise figures were used. Throughout this paper, a disc 
model with an accretion rate of 10 _7 Mq yr _1 is used, 
but in Table El we also present results of time-dependent 
chemical models for a lower accretion rate, M = 10 -9 M 
yr -1 . 

To characterise differences between model predictions, 
we divide all species in three groups by the maximum 
ratio of their theoretical column densities (less than 20, 
less than 200, and the rest). Given the wide variety of 
conditions and assumptions made in the various models, 
it is natural to expect that differences in calculated column 
densities should be significant. However, in reality, column 
densities for many species are close to each other in the 
different studies. For example, the scatter in HCN column 
densities does not exceed an order of magnitude. Similarly, 
ammonia densities in all but one case are a few times 
10 13 cm -2 . Note that we omitted the anomalously high 
NH3 column density computed by Aikawa & Herbst ( 1999 ) 
as this was based on an unrealistic rate for the H^" + N 
reaction (Aikawa & Herbst 20OTJ. 

Somewhat higher maximum-to-minimum ratios for 
species in the second group are mainly caused by their low 
abundances in the Willacy & Langcr ( 2000 ) model. In the 
last column of Table we give relative variations of their 
column densities both with and without (in parentheses) 
Willacy & Langer ( 2000) data. If the latter are not taken 
into account, the ratio of maximum to minimum column 
density does not exceed 50 for all these molecules. This 
similarity of column densities and a lack of correlation be- 
tween them and the column density of molecular hydro- 
gen can be understood as a further manifestation of the 
"warm molecular layer" (Aikawa et al. I2002[) where nearly 
all admixture gas-phase molecules are concentrated. On 
the other hand, molecular hydrogen is concentrated in the 
midplane, where all other molecules are frozen out and 
thus do not contribute to column densities. 

The reason why column densities in the calculations 
of Willacy & Langer are lower in comparison with other 
studies is probably related to the disc model they adopted. 
In the model of Chiang & Goldreich l|1997|) the disc is as- 
sumed to consist of only two layers, namely, the dark dense 
midplane and the less dense surface layer subject to harsh 
UV radiation. In the cold midplane molecules are mainly 
locked in icy mantles while in the surface layer they are 
easily dissociated or ionised by UV photons. In contrast 
to the disc model of DAlessio et al. 1)19990 . there is no 
region similar to the shielded intermediate ( "warm molec- 
ular") layer, which is the most favourable disc domain for 
many molecules to have their maximum gas-phase concen- 
trations. 



The most striking finding is variations in theoreti- 
cal CO column densities. As carbon monoxide, like other 
molecules, has its highest abundance in the intermediate 
layer, we might expect that its abundance is relatively in- 
dependent of the total gas column density. Thus, the ori- 
gin of the large discrepancies in CO abundances should lie 
in physical differences between various models. Looking at 
Table 03 one is tempted to assume that the CO abundance 
is low in models with surface chemistry (this paper and 
Willacy & Langer 2000 i or in models with X-ray induced 
chemistry (this paper and Aikawa & Herbst IT9991 1200 
In the former case the CO molecule might be transformed 
to CO2 and H2CO, in the latter case it might be destroyed 
by abundant He + . 

To check if this is the case, we computed the vertical 
distribution of CO for R = 370 AU without surface reac- 
tions and with Cx=0. In both cases we failed to reproduce 
the high CO column densities obtained by Aikawa et al. 
( 20021 and van Zadelhoff et al. (200$. As we noted above, 
the low CO column density in the Willacy & Langer model 
can partly be caused by the disc model they adopted. 
A possible explanation for the low CO abundance in our 
model is the fact that we do not take into account self- 
and mutual-shielding of molecular hydrogen and carbon 
monoxide, contrary to Aikawa et al. I|2002J) and Zadelhoff 
et al. <j2T)T73)l . 

In addition, we include in Table observationally in- 
ferred column densities for DM Tau (single-dish measure- 
ments) and LkCal5 (interferometric data). Note that DM 
Tau values correspond to column densities averaged over 
the entire disc (~ 800 AU), while values estimated from 
LkCal5 single-dish observations are column densities av- 
eraged over a ~ 100 AU disc. For species in the first and 
second groups there is a reasonable agreement between 
our predictions and observational data. 

6. Discussion 

Recent theoretical studies (e.g. Aikawa & Herbst 1999 
I21M1 Aikawa et al. I2TKH van Zadelhoff et al. 2003(1 re- 
vealed that the distribution of gas-phase molecular abun- 
dances within a steady-state accretion disc in the absence 
of mixing processes has a three-layer vertical structure. In 
the dense and cold midplane all species, except the most 
volatile ones, are frozen out on grain surfaces while in the 
disc atmosphere (surface layer) they are easily destroyed 
by high-energy stellar radiation. Therefore, in these disc 
domains gas-phase abundances of polyatomic species are 
low. On the contrary, in the intermediate layer, which is 
shielded from UV photons but still warm enough to al- 
low effective desorption, most molecules remain in the gas 
phase and drive a complex chemistry resulting in a rich 
molecular composition. As we noted above, this simple 
picture may not be appropriate when diffusion processes 
are taken into account (Ilgner et al. 2003). 

In our study we focused on the analysis of chemical 
processes relevant to the evolution of disc fractional ioni- 
sation based on the reduction approach. The size of a re- 
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Table 9. The observed and calculated column densities in cm 2 , r = 370 AU, t — 10 6 yrs. 



Species 


I 


II 


III 


IV 


V 


This 


paper 




Observed 




Max. 














M, M Q yr" 1 


DM Tau 


LkCal5* 


LkCal5** 


difference 














10" 9 


10- 7 


(SD) 


(IM) 


(SD) 




H 2 


9(22) 


9(21) 


1.5(24) 


1.3(22) 


1.3(23) 


4.8(21) 


1.4(24) 


3.8(21) 








H 2 


_ 







2.7(14) 


8.0(13) 


2.1(14) 


1.6(14) 











~ 3 


HCN 


2(12) 


2(12) 


2.1(12) 


1.4(12) 


2.5(12) 


1.3(13) 


6.7(12) 


2.1(12) 


0.02 - 1.2(15) 


7.8(13) 


~ 10 


NH 3 




1(13) 





1.4(13) 




4.4(13) 


1.4(14) 









~ 15 


C 2 H 


2(12) 


8(12) 


6.2(12) 


8.9(11) 


2.0(13) 


7.0(12) 


1.1(13) 


4.2(13) 








~ 20 (10) 


HNC 


1(12) 




2.0(12) 


3.0(11) 





6.9(12) 


8.0(12) 


9.1(11) 


< 5.4(12) 





~ 30 (8) 


CN 


1(13) 


1(13) 


3.8(12) 


5.3(11) 


3.0(13) 


6.4(12) 


4.6(12) 


9.5 - 12(12) 


9.1 - 25(13) 


6.3(14) 


~ 60 (8) 


N 2 H+ 






1.9(12) 


8.1(11) 




6.0(12) 


8.1(13) 


< 7.6(11) 


< 5.7(12) 


< 5.9(13) 


~ 100 (40) 


HCO+ 


3(11) 


5(11) 


9.0(12) 


8.8(10) 


9.0(12) 


8.4(11) 


2.0(11) 


4.6 - 28(11) 


1.5(13) 


1.4(13) 


~ 100 (45) 


cs 


3(H) 


1(12) 


4.9(11) 


1.0(11) 


1.2(12) 


5.6(12) 


1.1(13) 


6.5 - 13(11) 


1.9-2.1(13) 


2.2(14) 


~ 110 (35) 


CO 


2(15) 


8(15) 


1.1(18) 


7.1(16) 


1.2(18) 


1.6(16) 


2.5(16) 


5.7(16) 


1.6(18) 


9.0(17) 


~ 600 


H 2 CO 


7(11) 


1(12) 


2.9(12) 


8.7(10) 




8.0(13) 


8.3(13) 


7.6-19(11) 




3.0-22(13) 


~ 10 3 


CH 3 OH 






6.4(08) 


7.7(11) 




4.0(10) 


2.9(10) 




7.3 - 18(14) 


< 9.4(14) 


< 10 3 



I - Aikawa & Herbst (1999), Fig. 7, R = 400 AU, t = 9.5 • 10 5 yrs, high £ case 

II - Aikawa & Herbst (2001), Fig. 6, R = 400 AU, t = 9.5 • 10 5 yrs 

III - Aikawa et al. llWfll . Table 1, M = 10" 7 M e yr -1 , R = 370 AU, t = 9.5 • 10 5 yrs 

IV - Willacy & Langer (2000), Table 4, interpolated to R = 400 AU, t = 10 6 yrs 

V - van Zadelhoff et al. $Mty . Fig. 5, R = 400 AU 
IM - interferometric observations, beam size is ~ 4" 
SD - single-dish observations, beam size is ~ 20" 

* - do not necessarily correspond to r = 370 AU, 
** - disc size is ~ 100 AU 



duced chemical network, accurately reproducing the elec- 
tron concentration in a given disc region, can be under- 
stood as a quantitative criterion for the complexity of ion- 
isation chemistry. In Figure [7] we show the distribution of 
the number of species in the reduced networks over the 
disc. This distribution demonstrates a layered pattern. In 
most parts of the disc the chemistry of ionisation is sim- 
ple either due to the lack of ionising factors, low temper- 
ature and high density (midplane) or due to the presence 
of ionising radiation (surface layer). In these regions, it 
is sufficient to keep about 10-30 species and a few tens 
of reactions in reduced networks. Between the midplane 
and the surface layer the intermediate layer is located, 
where the evolution of ionisation degree is more compli- 
cated from the chemical point of view (especially at early 
times). Here, one has to retain more than fifty species 
and a comparable number of reactions in the reduced net- 
works. 




1 10 100 

R, AU 



Fig. 7. The distribution of the number of species in the 
reduced networks governing fractional ionisation over the 
disc. The highest and smallest number of species is 112 
and 5, respectively. 



This layered structure is directly related to the size and 
location of the so-called "dead zone" within a protoplane- 
tary disc, which is the region where the ionisation degree 
is so low that the magnetic field is not coupled to the 
gas. The MHD turbulence cannot develop in this region. 
Thus, the MRI is not operative, which means the lack 
of an effective angular momentum transport if no other 
transport mechanism can be identified. The "dead zone" 
in accretion discs has been widely investigated (e.g. Igea 
& Glassgold IT3531 Fromang et al.EEEl Sano et al. TMKty . 



An important quantity which characterises the cou- 
pling between the matter and the magnetic field is the 
magnetic Reynolds number, which can be defined as 
(Fromang et al. 127)172^1 



Re 



M 



c s H 

•j 

V 



(7) 



where 77 is the Ohmic resistivity (Blaes & Balbus fl994|l 
234 ■ T°- 5 

rj= . (8) 
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Gas-grain chemistry, 
t = 10" yrs 




3.8e + 013 



10 



100 



7.9e + 000 



R, AU 



Fig. 8. The magnetic Reynolds numbers computed for our 
model and an evolutionary time of 10 6 yrs. The solid line 
designates Reffi = 100. The highest and lowest values are 
equal to about 10 13 and 10, respectively. 



The other quantities are the sound speed c s , and the thick- 
ness of the disc H . The MHD turbulence can be main- 
tained only if this quantity exceeds a certain critical value 
which depends on the field geometry and other factors. 
Following Fromang et al., we consider two cases of this lim- 
iting value, namely, Reffl = 100 and i2e •?* = 1 and define 
the "dead zone" as a disc domain where Reu < R^m"- 

The calculated magnetic Reynolds numbers for our 
model are shown in Fig. [5] The lowest magnetic Reynolds 
number is R&m ~ 10, which implies that under certain 
circumstances within the framework of our model the 
"dead zone" is entirely absent. With another critical value, 
flej 1 ' = 100, the "dead zone" occupies the following disc 
region: 2 AU < r < 20 AU, z ~ 0.06 • r 1 - 27 (solid line). 
This result is roughly consistent with other recent studies. 
For example, Fromang et al. found that the "dead zone" 
can extend over 0.7 AU < r < 100 AU, if the viscos- 
ity parameter is equal to 10 -2 and the accretion rate is 
10 Mg yr — x , whereas Igea & Glassgold (1999) obtained 
a "dead zone" of somewhat smaller size, r < 5 AU. 

Apart from defining the location of the "dead zone" 
in accretion discs, the fractional ionisation determines 
whether non-ideal MHD effects are important for the evo- 
lution of protoplanetary discs. When the electron abun- 
dance is computed for the solution of non-ideal MHD 
equations, an ionisation equilibrium is often assumed. An 
expression similar to Eq. (JJJ is then used to calculate the 
equilibrium fractional ionisation, with different estimates 
for the gas-phase recombination coefficient (3 (Blaes & 
BalbusESl Ganmiielinni Regos HWfl Rey es-Ruiz I2"0lfl1 
Fromang et al- EHoU Fleming & Stone EDD5|l . 

In Fig. we compare the fractional ionisation x c (eq), 
computed from Eq. Q with the often used expression for 
recombination coefficient (3 — 8.7 • 10~ 6 T -1 / 2 (Glassgold, 
Lucas, & Omont fl986[) . to the fractional ionisation com- 
puted with the full chemical network x c (full), for some 
representative points in the midplane and in the inter- 
mediate layer. It is not surprising that the ratio of these 
two quantities exceeds 10 at low ionisation degrees. Here 



8" io° 




-14 -13 -12 -11 -10 -9 
log x.(full) 



-8 -7 



Fig. 9. Comparison of the equilibrium and time- 
dependent fractional ionisations at t = 10 6 years in the 
midplane and in the intermediate layer. Symbols corre- 
spond to different radii. 



the electron density is small and ion-grain interactions 
should be taken into account. At moderate ionisation de- 
grees (10~ 10 — 10~ 8 ), the equilibrium value may differ from 
the "exact" value by a factor of a few at the end of the 
computation. At earlier times, the discrepancy is higher. 

However, in the midplane the equilibrium state is 
reached very rapidly, with no appreciable changes in x c 
after ~ 1000 years of evolution. If one is not interested in 
time scales less than 1000 years and can afford a factor 
of a few uncertainty in the fractional ionisation, deep in 
the disc interior the equilibrium x c is sufficient, provided 
grain charging processes are taken into account. 

The intermediate layer shows a very different be- 
haviour. As we mentioned previously, the equilibrium 
state is not reached there or is reached very late. The ini- 
tial fall-off of x e , caused by HCO + , Hjj" and metal ion re- 
combination, is followed by a step-like increase. The dura- 
tion of each phase and x c evolutionary changes, illustrated 
in Figure^] vary significantly in different parts of the layer. 
Surprisingly, the equilibrium x e from Eq. is very close 
to the computed fractional ionisation at t — 10 6 years in 
this point (1.5 • 10~ 9 and 2 • 10~ 9 , respectively). This is, 
of course, a coincidence, as the equilibrium is not reached 
there. We emphasise that even if z c (full) is equilibrated 
eventually in the intermediate layer and is close to x (eq), 
the latter value does not reflect the ionisation state of the 
medium for most of the computational time. 

In the surface layer the fractional ionisation is con- 
trolled not only by X-rays, but also by UV-radiation. 
The equation is still applicable here but with an im- 
portant change. In the surface layer x e is determined 
mostly by ionised carbon (Table Its recombination 
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coefficient is much lower than the value quoted above, 
thus, /3 = 8.7 • lO^T- 1 / 2 should not be used in this 
case. If one substitutes in Eq. Q the correct /3 for C + 
(1.4 • lCr 13 T-°- 61 ), the resultant equilibrium few 
times smaller than the "exact" value, mainly because pho- 
toionisation is not taken into account. Time needed to 
reach the equilibrium x e in the surface layer does not ex- 
ceed 100 years. It is questionable whether precise informa- 
tion about such a high ionisation (x c ~ 10~ 4 ) is required 
for MHD computation. 

Overall, we conclude that the equilibrium approach is 
appropriate most of the time and in most parts of the 
studied disc, but it should be applied with care. This is es- 
pecially important in less dense discs that are more trans- 
parent to X-rays and UV radiation. In this paper, we used 
a relatively massive disc model with a narrow intermediate 
layer. In models with lower accretion rates the intermedi- 
ate layer would occupy a greater volume of the disc. 

Even though the reduced networks presented in this 
paper provide much better accuracy than the equilibrium 
approach, while having quite modest sizes, their usefulness 
for dynamical models is limited by the fact that they are 
computed in a steady-state disc. For example, in the point 
Ml the reduced network consists of three species, NO, 
Mg, and Mg + . This network is of no use in a dynamical 
model because there will be no NO in the gas flowing 
into the inner part of the midplane, as reduced networks 
further out in the midplane do not contain this molecule. 
To make the reduction valuable for MHD modelling we 
need to consider it in a dynamically evolving medium. 

An alternative approach is to merge all the networks 
described in this paper into a single network that com- 
prises about 120 species and a comparable number of re- 
actions. This network is out of the scope of the present 
paper and will be described in further publications. 

7. Conclusion 

Chemical processes, responsible for the ionisation struc- 
ture of a protoplanetary disc with a central star, are anal- 
ysed by means of reduced networks that reproduce the 
ionisation fraction within a factor of 2. These networks 
are available from the authors upon request. Because of 
the wide range of physical conditions met in a typical 
disc, there is a corresponding diversity in the chemical 
reactions that control the fractional ionisation in differ- 
ent parts of the disc. Generally, it can be divided into 
three layers. In the midplane the ionisation is provided 
by cosmic rays and radioactive elements only. Above the 
midplane, the intermediate layer is located where the ion- 
isation is dominated by X-rays. In the surface layer UV 
photons are the main ionising factor. In each of these lay- 
ers we analyse several representative points and construct 
reduced chemical networks that are needed to reproduce 
the fractional ionisation as a function of time during 10 6 
years of evolution within a factor of 2 uncertainty. In the 
midplane the chemistry, which determines the fractional 
ionisation, is very simple. Reduced networks typically in- 



clude no more than 10 species and a similar number of 
reactions. The same holds true for the surface layer. The 
intermediate layer is the most complicated region from the 
chemical point of view. To compute the fractional ionisa- 
tion with the targeted accuracy, in some regions one has 
to take into account carbon chains containing up to 6 car- 
bon atoms, which leads to reduced networks with over a 
hundred species and reactions. 
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